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The  bispectrum  is  the  natural  third-order  generalisation  of  the  power  spectrum.  It 
provides  information  about  correlations  between  different  Fourier  components  of  a 
signal  or  image,  and  about  the  statistics  of  Fourier  phase.  A number  of  numerical 
and  experimental  studies  of  the  bispectra  of  chaotic  systems  have  been  published. 

In  this  paper  we  present  the  first  analytical  calculations  of  the  bispectra  of  chaotic 
dynamical  systems.  First,  for  a generalisation  of  the  classical  sawtooth  or  Renyi 
map,  we  calculate  the  bispectrum  using  symbolic  dynamics.  Also,  for  intermittent 
systems,  we  calculate  the  bispectrum  using  the  relationship  between  these  systems 
and  renewal  processes.  We  review  the  results  of  these  calculations,  drawing  some 
conclusions  about  the  characteristic  features  of  the  bispectra  of  chaotic  systems, 
and  compare  them  with  the  features  of  some  financial  time  series. 


Key  words:  Chaos,  bispectrum,  time  series,  intermittent,  finance. 

1 Introduction:  power  spectrum  and  phase  in  time-series  analysis 

Given  an  experimental  time  series  xt,  what  can  we  say  about  the  nature  of  the 
system  that  produced  it?  This  is  a question  commonly  encountered  by  experimental 
scientists,  and  also  by  those  who  deal  with  financial  time  series  such  as  stock  prices 
or  currency  exchange  rates.  Often  the  first  attack  at  the  question  is  to  calculate  the 
power  spectrum  P{uf).  The  power  spectrum  reveals  how  the  variance,  or  power,  in 
the  time  series  is  shared  between  the  different  frequencies  uj. 

Because  the  power  spectrum  contains  no  phase  information,  it  is  not  possible  to 
reconstruct  the  original  signal  xt  from  it,  and  there  are  some  questions  that  can  not 
be  answered  by  inspecting  only  the  power  spectrum.  One  of  these  is  the  question 
of  whether  an  apparently  random  signal  originates  from  chaos  in  a low-dimensional 
deterministic  dynamical  system  or  from  a stochastic  process.1,2 
The  Fourier  transform 

= J2Xtetut  (!) 

t 

contains  the  same  information  as  the  original  signal,  and  x can  be  written  in  terms 
of  a modulus  and  a phase: 

*M  = |*(w)|e*<"\  (2) 

where  the  modulus  |x(w)|  determines  the  power  spectrum.  To  go  beyond  the  power 
spectrum  while  still  using  the  concepts  of  Fourier  analysis,  one  should  therefore 
consider  the  statistics  of  the  Fourier  phase  0(o>).  In  fact,  the  phase  is  more  impor- 
tant in  determining  the  appearance  of  an  image  than  the  Fourier  modulus  or  power 
spectrum.3  All  information  about  the  location  of  features  in  an  image  or  the  time 
of  occurrence  of  features  in  a signal  is  contained  in  the  phase. 
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It  is  not  straightforward  to  study  the  statistics  of  phase  in  a time  series.  One 
advantage  of  the  power  spectrum  is  that  it  is  a well-defined  statistical  average, 
being  the  Fourier  transform  of  the  autocorrelation  function.  However,  stationary 
stochastic  processes  (including  those  defined  by  a chaotic  dynamical  system  with 
its  invariant  measure)  do  not  possess  a well-defined  average  phase  cj>(u>)  or  e 
for  any  frequency  u. 

We  can  understand  why  the  Fourier  phase  factor  el<t> ^ does  not  have  an  average 
value  for  a stationary  process  by  noting  that  translating  the  original  time  series  by 
a delay  s (that  is,  replacing  xt  by  X(t_s))  is  equivalent  to  multiplying  each  Fourier 
phase  factor  by  eVjJS . This  is  just  a restatement  of  the  ‘shift  theorem’  of  Fourier 
analysis.  Since  any  time-average  of  a stationary  process  must  be  invariant  under  a 
time-shift,  this  explains  why  the  phase  factor  does  not  have  a well-defined  average. 

This  argument  also  provides  a hint  of  how  we  might  construct  well-defined 
average  quantities  that  contain  phase  information.  For  any  two  frequencies  oq  and 
u> 2,  the  phase  product 

R(u>  i,wa)  =ew"l)e*("a)e_*{"1'H*)  (3) 

is  invariant  under  a translation.  A non-zero  average  of  this  quantity  implies  that  the 
phases  of  the  Fourier  components  at  frequencies  oq,  w2  and  oq  + w2  are  correlated. 

In  fact,  the  mean  value  of  R is  the  phase  of  the  bispectrum, 4,5  evaluated  at  the 
frequencies  (nq,uq).  Just  as  the  power  spectrum  P{ u)  is  the  Fourier  transform 
of  the  autocorrelation  function,  the  bispectrum  P(uq , oq)  is  the  double  Fourier 
transform  of  the  bi-correlation  function 


Cx{s,  t)  = ((xs  - x ) (xt  - a;)  (x0  - x)) , (4) 

where  the  angled  brackets  ()  denote  an  expectation  value. 

The  bispectrum  can  also  be  estimated  directly  from  the  Fourier  transform: 

P{loi,u2)  = Jim  ~ { xt{uii)xt{u}2)x*t{uji  + w2)> , (5) 

where  the  subscript  T implies  that  the  sum  in  equation  (1)  is  over  the  first  T 
samples  in  the  time  series,  and  =1=  denotes  complex  conjugation. 

Higher-order  generalisations  of  the  power  spectrum,  known  as  polyspectra,  can 
also  be  defined4’5.  The  bispectrum  has  been  applied  in  time-series  analysis  to  con- 
struct tests6, 7,8for  nonlinearity  and  for  normality  (meaning  the  presence  of  a normal 
distribution).  These  tests  depend  upon  the  observations  that  for  normal  systems, 
the  bispectrum  is  zero,  and  for  linear  systems  (that  is,  processes  constructed  by  a 
linear  operator  acting  upon  an  i.i.d.  process)  the  bicoherence 


&(uq,w2)  = 


|P(uq,w2)|2 


P(uq)P(u;2)P(uq  + w2) 


(6) 


is  constant. 

The  bispectrum  has  also  been  used  to  study  phase  correlations  and  the  in- 
teractions between  Fourier  components  in  many  experimental  systems,  including 
electroencephalographic  signals  from  the  human  brain,9  and  the  oscillations  of  the 
earth.10  A number  of  authors8,11,12,13  have  studied  the  bispectra  of  time  series 
from  chaotic  dynamical  systems,  either  by  numerical  simulation  or  experimental 
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Figure  1.  Symmetries  of  the  bicorrelation  function.  For  example,  if  the  bicorrelation  is  given  by 
f(s,t ) in  the  region  A,  then  the  bicorrelation  in  region  C is  given  by  f(t  — s,s ). 


measurement.  Subba  Rao8  identified  ‘ridges’  in  the  bicoherence  as  a characteristic 
of  chaotic  systems,  and  other  authors  studying  electrical  circuits  have  noted  vari- 
ous changes  in  the  bispectral  characteristics  of  time  series  as  the  systems  generating 
them  make  the  transition  to  chaos. 

To  date,  there  are  no  published  examples  of  chaotic  dynamical  systems  where 
the  bispectrum  can  be  calculated  exactly.  In  this  paper  we  present  an  example  of 
such  a calculation,  for  a simple  discrete-time  system  which  is  a generalisation  of  the 
well-known  saw-tooth  or  Renyi  map.  We  also  show  how  the  relationship  between 
intermittent  chaotic  systems  and  renewal  processes  can  be  used  to  give  asymptotic 
low-frequency  expressions  for  the  bispectra  of  intermittent  systems.  We  identify 
some  common  features  of  these  two  families  which  may  be  bispectral  character- 
istics of  chaos.  Finally,  we  compare  some  financial  time  series  with  the  chaotic 
examples.  Before  presenting  the  calculations,  we  briefly  review  some  properties  of 
the  bicorrelation  function  and  bispectrum. 


2 Symmetries  of  the  bicorrelation  function  and  calculation  of  the 
bispectrum 

This  section  summarises  some  properties  of  the  bicorrelation  function  and  bispec- 
trum that  are  useful  in  deriving  the  form  of  the  bispectrum.  From  its  definition  in 
equation  (4),  one  can  see  that,  for  a stationary  process,  the  bicorrelation  function 
Cx(s,t)  has  the  symmetries 

Cx(s,t)  = Cx(t,s)  Cx(~s,t)  - Cx(s,s + t).  (7) 

By  composing  these  symmetry  operations,  one  can  reach  any  part  of  the  (s,  t ) plane 
from  the  region  0 < s < t where  the  bicorrelation  is  most  easily  calculated.  These 
relationships  are  summarised  in  figure  1. 

The  bispectrum  P(uq,o; 2)  is  defined  as  the  double  Fourier  transform  of  the 
bicorrelation  function,  the  Fourier  integral  being  evaluated  over  the  entire  ( s,t ) 
plane.  Let  Pa(^i,  W2)  be  the  integral  over  only  region  A in  figure  1.  If  the  value  of 
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f(x) 


f(x) 


Figure  2.  The  two  types  of  dynamical  system  considered  in  this  paper:  the  generalised  sawtooth 
map  (left)  with  b = 0.3,  and  the  intermittent  dynamical  system  described  by  equation  (28)  with 
fi  = 1,  z = 2.7. 


the  bicorrelation  in  region  A is  given  by  Cx(s,t ) = 0(s,t),  then 

r-OO  rOO 

Pa(oj1,uj2)=  ds  / dt  /?(M)ei(sWl+ta,2). 
Jo  j s 

Similarly,  the  integral  over  region  B is 

/0  r co 

ds  / dt  Cx(s,t)e^su,1+t^l 

-co  J 0 


(8) 

(9) 


Using  the  symmetries  in  equation  (7),  and  changing  the  variables  of  integration, 
one  can  show  that 


PB{ui,u)2)  = PA{-{vi +u>2),lj2).  (10) 

Similar  relationships  can  be  derived  between  PA  and  the  Fourier  integrals  over  the 
other  regions  in  figure  1.  The  bispectrum  can  then  be  calculated  by  summing  all 
the  integrals.  The  result  is 

P(u>i,oj2)  = PA(uii,  u)2)  + Pa{~{u\  + uj2),u;2)  + PA(u)i,  —(u>i  + u2)),  (11) 

where  the  hat  operation  ? denotes  symmetrisation: 

f{x,y)  = f{x,y)  + f(y,x).  (12) 


3 The  generalised  sawtooth  map 

3.1  The  dynamical  system 

Our  first  dynamical  system  is  defined  by  the  mapping  of  the  interval  (0, 1)  onto 
itself 


x < b 
x > b 


f{x)  = ( x-6 

1 i-i> 

as  pictured  in  Figure  2.  The  dynamical  system  xt+i 


(13) 

f(xt ) has  a uniform  invariant 
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distribution.  It  is  chaotic,  with  Lyapounov  exponent 

A = -blnb  - (1  - b)  ln(l  - 6)  (14) 

easily  calculated  by  standard  methods.14 

This  dynamical  system  is  a generalisation  of  the  sawtooth  or  Renyi  map  which 
has  been  studied  by  many  authors.14  When  b = 1/2,  the  two  systems  are  identical. 
This  generalisation  makes  the  study  of  the  system  more  difficult,  but  it  is  necessary 
in  this  case  because  when  b = 1/2,  the  system’s  bispectrum  is  zero.  This  is  an 
effect  of  the  symmetry  of  the  function  /.  Sakai  and  Tokumaru15  have  calculated 
the  autocorrelation  function  of  the  6 = 1/2  system  using  a relationship  between 
this  system  and  an  autoregressive  process. 

3.2  Symbolic  dynamics 

The  method  we  use  to  study  this  system  is  symbolic  dynamics  14 , where  a sequence 
of  symbols  from  a discrete  and  finite  alphabet  represents  a trajectory  of  a dynamical 
system.  The  representation  for  the  generalised  sawtooth  map  is  not  difficult  to 
derive.  Suppose  that  xq  and  x\  are  two  successive  points  in  the  trajectory  of  the 
system,  so  that  x\  = f(x 0).  Then  we  have  either 

xq  = bxi  (15) 

or,  if  x0  > 6, 

xq  = b + (1  - b)xi.  (16) 

Putting  z\  = 0 in  the  first  case,  and  Z\  = 1 in  the  second,  we  can  write 

Xq  = z\b  + x\e\,  (17) 

where  e*  = (6  - 2 z$  + Zi).  Continuing  along  the  same  lines,  we  obtain 
xq  - zib  + ei(z2b  + e2(z3b  + e3x3)) 

n 

= 6eie2---eTl_ixn  + 6]r]eie2---efc_i;Zfc.  (18) 

k= 1 

It  is  easy  to  show  that,  under  the  (uniform)  invariant  measure,  the  Zi  are  indepen- 
dent random  variables,  taking  the  values  0 or  1 with  probabilities  b and  (1  — 6) 
respectively.  We  can  therefore  use  this  representation  to  calculate  various  averages. 
For  example,  in  the  next  section  we  calculate  the  autocorrelation  function  and  the 
power  spectrum. 

3.3  The  power  spectrum 

Multiplying  equation  (18)  by  xn,  we  obtain 


n 

x0xn  = beie2  ■ ■■en-.ix1n  + 6^eie2  ■ • -ek^izkxn. 

k= 1 


(19) 
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We  can  now  take  the  expectation  value,  using  the  fact  that  random  variables  with 
different  subscripts  are  independent,  and  the  expectation  values  (zi)  = (1  — 6), 
(zi)  = (1  - 6),  (xi)  = 1/2,  (xj)  = 1/3  and 

(ei)  = 62(1  - 6)2  = B,  (20) 

(e2)  = 63  + (1  - 6) 3 = (7,  (21) 

(22) 

After  summing  a geometric  series,  we  obtain 

1 Bn 

( XoXn ) — ^ + -J2'  (23) 


This  equation  is  valid  for  n > 0.  The  autocorrelation  function  is  therefore 


Cx(t)  = ( x0xt ) - (x20)  = (24) 

The  power  spectrum  is  the  Fourier  transform  of  this  function  (using  the  convention 
of  equation  (1)): 


P(uj)  = 


1 + B2 

12(1  + B2  —2B cos#)  ’ 


(25) 


3-4  The  bispectrum 

The  same  methods  can  be  used  to  calculate  the  three-point  average  {x03. :sxt)  when 
0 < s < t.  We  can  then  calculate  the  bicorrelation  function.  The  details  of  this 
calculation  occupy  far  too  much  space  to  be  recorded  here.  The  result  is,  again  for 
0 < s < t, 


Ct(s,  t)  = {Bl  - i?5^-5)  , (26) 

where  B and  C are  defined  in  equations  (20)  and  (21).  The  bispectrum  can  now 
be  calculated  using  the  method  described  in  the  last  section.  The  result  is  given  by 
equation  (11),  with 


Pa{x,v) 


(26  — 1)(6  — l)be~ix 
12  (1  - Be~ix ) (1  - Ce~ix)  (1  - Be^) 


(27) 


Contour  diagrams  of  the  bispectrum  and  bicorrelation  function  for  a typical 
value  of  6 are  shown  in  figure  3.  The  ‘ridges’  found  by  Subba  Rao8  are  not  evident 
for  this  system  for  any  value  of  6 that  we  have  examined.  The  most  noticeable 
consistent  feature  is  that  the  energy  in  the  bispectrum  is  spread  over  a broad 
ranges  of  frequencies  rather  than  being  concentrated  into  narrow  bands.  This  may 
be  a characteristic  of  fully-developed  chaos. 
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Figure  3.  The  bispectrum  for  the  sawtooth  map  with  b = 0.27.  The  left  graph  shows  a surface 
map  of  the  modulus  of  the  bispectrum,  evaluated  over  the  region  where  ui i and  U2  are  positive. 
The  right  graph  shows  a direction  field  displaying  the  complex  phase  of  the  bispectrum  overlaid 
with  a contour  map,  evaluated  over  both  positive  and  negative  frequencies. 


4 Intermittent  dynamical  systems 


In  this  section,  we  study  the  bispectra  of  a class  of  intermittent  dynamical  systems. 
For  definiteness,  we  take  as  our  model  the  system  defined  by  the  mapping 

f{x)  — (x  4-  fj,xz)  mod  1,  (28) 


shown  in  figure  2 for  /r  = 1,  2 = 1.7.  The  analysis  given  here  depends  only  on  the 
form  of  f(x)  close  to  the  intermittent  point  at  x = 0,  and  so  applies  to  a broad 
class  of  intermittent  dynamical  systems,14  as  described  by  Ben-Mizrachi  et  al.16 
The  defining  characteristic  of  these  intermittent  systems  is  that  they  show  brief 
periods  of  random  behaviour  between  long  periods  when  the  system  is  ‘stagnant’, 
remaining  close  to  the  centre  of  intermittency  (in  the  case  of  the  model  system 
above,  close  to  the  point  x = 0). 

We  will  be  studying  the  long-time  behaviour,  or  equivalently  the  low-frequency 
behaviour,  of  these  dynamical  systems.  For  this  purpose,  following  Ben-Mizrachi  et 
al,16  we  replace  the  dynamical  system’s  time  series  xt  by  a series  of  delta-functions 
located  at  the  times  of  escape  from  periods  of  stagnation.  The  dynamical  system 
is  thus  replaced  by  a renewal  process  with  waiting-time  distribution  defined  by  the 
distribution  of  waiting  times  between  escape  events.  As  Ben-Mizrachi  et  al 16  show, 
this  distribution  can  be  approximated  for  large  r by 


Pt 


{a  - 1 )t~a 
0 


t > 1 
t < 1, 


(29) 


where  a = z/(z  — 1).  Ben-Mizrachi  et  al 16  also  show  that  the  function  c(t)  giving 
the  probability  of  an  escape  event  at  time  t given  the  occurrence  of  such  an  event 
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at  time  0 is  given  by 


c(s)  = 


1 

1 -Pr(s)’ 


(30) 


where  the  bar  . denotes  the  Laplace  transform.  The  Laplace  transform  of  pT  is  an 
incomplete  Gamma  function,1'  whose  small-s  behaviour  is  given  by 


1 — pT(s)  ~ s“_1r(2  — a) 


(31) 


for  1 < a < 2.  For  the  sake  of  brevity,  we  consider  only  this  case  (corresponding 
to  z > 2)  in  this  paper.  Other  cases  will  be  covered  in  a longer  paper  currently  in 
preparation. 

For  the  case  1 < a < 2,  for  small  s we  have  c(s)  ~ T(2  — a)s1-“,  and  the 
low-frequency  behaviour  of  the  power  spectrum16  is  given  by 

P(w)  = I c(iu)  + c(-iu) I2  ~ ^2a_2  ■ (32) 

The  bispectrum  can  also  be  calculated  from  c(s).  First,  note  that  for  large  t\ 
and  (f2  — h ) and  for  0 < t\  <t2, 

(xtlx0)  ~ c(oo)c(ti).  (33) 

and 


(xtlxt2x0)  ~ c(oo)c(ti)c(i2  - ti).  (34) 

The  bicorrelation  Cx(ti,t2)  is  therefore  given  by 

Cx(ti,t2)  = c(oo)c(fi)c(f2  - ti)  -c(oo)2  (c(ti)  + c(f2)  + c(t2  - ti ))  + 2c(oo)3 . (35) 

After  applying  the  double  Laplace  transform  and  taking  the  limit  s ->•  0,  we  find 
that  the  first  term  dominates  and 


Cx(si,s2) 


c(  oo) 

T2(2  — a)s“_1  (s2  — si))“-1 ' 


(36) 


By  setting  s i = iu> i,  s2  — we  obtain  an  expression  for  the  bispectrum  integral 
Pa- 


Pa{u i,w2)  ~ 


c(oo) 

r2(2  - a)  (iu i (w2  - wi))“_1 


(37) 


The  bispectrum  is  then  given  by  equation  (11). 

Bispectra  calculated  numerically  from  time  series  generated  directly  by  the  dy- 
namical system  (28)  are  in  good  agreement  with  this  result,  although  a large  number 
of  data  points  is  required  for  the  average  to  converge  to  a smooth  function.  Figure 
4 shows  the  result  of  one  comparison. 

The  structure  of  these  low-frequency  bispectra  is  simple.  They  show  power-law 
singularities  at  the  lines  uq  = 0,  w2  = 0 and  uq  + cu2  = 0.  Like  the  bispectra  for 
the  sawtooth  map,  the  bispectra  for  intermittent  systems  (at  least  at  low  frequen- 
cies) show  no  isolated  peaks.  Intensity  extends  across  a wide  range  of  frequencies 
modulated  by  a power-law  envelope. 
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Figure  4.  The  theoretical  treatment  given  in  this  paper  predicts  that  the  function  P{ui i,0)  should 
diverge  as  as  un  — > 0.  This  graph  shows  the  theory  (solid  line)  and  values  of  a numeri- 

cally calculated  bispectrum  using  10000  data  points  of  a time  series  generated  by  the  dynamical 
system  of  equation  (28)  for  z = 2.1  (a  = 1.909). 


These  results  suggest  that  bispectral  measurements  may  be  useful  in  cases  where 
it  is  difficult  to  distinguish  power-law  noise  from  a dynamical  system16,2  from  noise 
from  a truly  random  source,  such  as  filtered  white  noise.  White  noise  retains  its 
random  phases  when  passed  through  a linear  filter,  and  so  has  zero  bispectrum. 

5 Bispectra  of  financial  time  series 

The  bispectra  of  economic  time  series  were  first  computed  by  Godfrey,18  who  was 
able  to  reject  the  null  hypothesis  of  linearity  for  a number  of  stock  price  time  series, 
but  did  not  discuss  other  aspects  of  the  bispectral  form.  For  comparison  with  the 
other  examples  presented  here,  we  calculated  the  bispectrum  of  the  log-increments 
yn  = \og(xn/xn~i)  of  the  Dow  Jones  Industrial  Average  stock  market  index,  eval- 
uated daily  over  the  period  January  1994  to  February  1999.  The  bispectrum  was 
estimated  from  the  bicorrelation  using  a Parzen  window,  as  in  the  work  of  Subba 
Rao  and  Gabr7.  The  result  (figure  5)  shows  broad  spectral  intensity,  as  for  the  two 
chaotic  examples  presented  here,  but  concentrated  more  at  higher  frequencies.  This 
shift  towards  higher  frequencies  is  probably  a result  of  taking  the  log-increments, 
which  is  effectively  a differentiation  process. 

6 Conclusion 

One  recognised  characteristic  of  chaotic  dynamical  systems  is  that  they  show  a 
continuous  power  spectrum,  with  intensity  distributed  across  a wide  range  of  fre- 
quencies. The  results  of  this  paper  show  that,  at  least  in  some  cases,  this  is  also 
true  of  their  bispectra.  This  provides  a way  of  distinguishing  a chaotic  time  series 
from  filtered  white  noise.  However,  other  classes  of  processes  share  this  charac- 
teristic of  broad-band  bispectral  intensity.  This  is  clear  because  the  bispectra  for 
the  intermittent  systems  in  this  paper  were  calculated  by  exploiting  the  similarity 
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Figure  5.  The  left-hand  graph  shows  the  modulus  of  the  bispectrum  for  the  log-increments  of  the 
Dow  Jones  industrial  average  stock  index  (values  shown  on  the  right). 


between  these  systems  and  renewal  processes.  A broad-band  bispectrum  can  not 
therefore  be  taken  as  an  unambiguous  sign  of  chaos.  A broad-band  bispectrum 
does  indicate  interactions  between  Fourier  modes  over  a wide  range  of  frequencies, 
and  is  therefore  a sign  that  the  system  studied  does  not  have  a simple  description 
in  the  Fourier  domain.  The  broad-band  character  of  the  bispectra  of  financial  time 
series  should  therefore  be  taken  as  an  indication  of  the  complexity  of  the  underlying 
processes. 

The  authors  would  like  to  thank  Abby  Evans,  Martin  Turner  and 
Ursula  Augsdorfer  for  help  with  the  production  of  this  paper. 
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